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Multiscale  Modeling  of  Stiffness,  Friction  and  Adhesion  in  Mechanical  Contacts 
FA9550-0910232,  Mark  O.  Robbins,  Johns  Hopkins  University 


The  mechanical  behavior  of  contacts  plays  a  critical  role  in  the  function  of  devices  from 
aircraft  engines  to  microelectromechanical  systems  (MEMS).  As  the  demands  for  performance 
from  these  devices  increase  or  their  size  shrinks,  traditional  continuum  models  of  contact 
behavior  become  inadequate.  Overcoming  these  limitations  will  require  improved  understanding 
of  interconnected  processes  on  a  wide  range  of  scales  from  atomic  interactions  at  interfaces  to 
surface  roughness  on  nanometer  to  micrometer  scales  and  to  deformation  and  loading  at 
macroscopic  scales. 

One  objective  of  the  proposed  research  was  to  develop  and  extend  multiscale  modeling 
methods  that  are  capable  of  bridging  between  atomistic  and  continuum  models.  The  first  part  of 
this  report  describes  a  Green’s  function  method  that  allowed  us  to  extend  the  scale  of  simulations 
from  ~1  OOnm  to  ~3pm. 

The  second  objective  was  to  use  these  multiscale  methods  to  improve  understanding  of 
aspects  of  contact  mechanics  that  are  important  to  Air  Force  applications.  For  example,  the 
mechanical  behavior  of  complex  assemblies  such  as  gas  turbines  is  strongly  influenced  by  the 
stiffness  of  mechanical  joints  and  the  dissipation  associated  with  partial  slip  in  these  joints.  At 
the  other  extreme  of  length  scales,  MEMS  devices  often  fail  due  to  undesired  adhesion  that  is 
affected  by  surface  roughness  and  interfacial  forces.  Energy  loss  due  to  friction  forces  is 
important  to  devices  at  all  scales. 

The  second  section  of  the  report  describes  studies  of  single  asperity  contacts  that  test  the 
applicability  and  limits  of  commonly  used  continuum  descriptions.  One  study  showed  that 
adsorbed  monolayers  present  on  most  surfaces  have  a  profound  impact  on  contact  properties.  A 
second  showed  that  thermal  fluctuations  lead  to  nonnal  pressures  that  are  comparable  to  the  ideal 
yield  stress.  This  makes  it  difficult  to  extend  continuum  definitions  of  contact  to  atomic  scales. 
The  third  study  looked  for  partial  slip  in  contacts  with  tips  up  to  2  pm  in  radius.  The  partial  slip 
predicted  by  continuum  models  for  these  scales  was  not  seen.  However  a  simple  scaling  law  was 
discovered  for  an  unexpectedly  strong  size  dependence  in  the  coefficient  of  friction. 

The  third  section  of  the  report  describes  studies  of  multi-asperity  contacts.  The  contact 
correlation  function  for  the  commonly  used  bearing  area  model  was  obtained  and  shown  to 
disagree  with  numerical  calculations.  Extensive  calculations  of  the  stiffness  of  multiasperity 
contacts  were  used  to  test  a  simple  scaling  relation  that  depends  on  a  few  experimental 
parameters.  Finally,  a  new  theory  of  adhesive  contact  between  rough  surfaces  was  developed  that 
can  capture  the  changes  in  contact  area  with  load,  adhesion  strength  and  other  parameters. 

The  report  concludes  with  a  list  of  papers  published  and  in  preparation,  and  a  list  of  talks  on 
funded  research. 

1.  ALGORITHMIC  DEVELOPMENT 

A  great  challenge  to  tribological  studies  is  that  the  response  is  determined  by  physics  on  a 
wide  range  of  length  scales.  Friction  and  adhesion  forces  result  from  interactions  between  atoms 
on  opposing  surfaces.  An  accurate  model  must  include  the  explicit  geometry  of  these  atoms.  On 
the  other  hand,  the  separation  between  atoms  is  controlled  by  the  elasto-plastic  response  below 
the  surface.  To  correctly  capture  subsurface  effects,  one  must  include  the  response  to  depths  D 


that  are  comparable  to  the  lateral  extent  along  the  surface  L.  Meeting  both  these  criteria  in  a 
single  calculation  presents  great  computational  challenges. 

Our  current  fully  atomistic  simulations  of  contact  are  limited  to  linear  dimensions  less  than 
lOOnm.  The  corresponding  number  of  atoms  is  up  to  50  million,  and  calculations  take  days  on  a 
hundred  processor  cores.  Part  of  the  proposal  was  to  search  for  algorithms  that  allowed  much 
larger  systems  to  be  studied.  The  key  idea  is  that  subsurface  deformations  vary  slowly  and  may 
be  modeled  with  more  efficient  coarse-grained  methods. 

The  method  that  we  found  most  effective  was  a  generalization  of  the  Green’s  function 
approach  of  Campana  and  Muser  (2006)  that  was  implemented  into  the  molecular  dynamics  code 
LAMMPS  by  Kong  et  al.  (2009).  Atoms  above  some  plane  are  treated  fully  atomistically.  The 
linear  elastic  response  of  atoms  below  this  plane  is  fully  captured  using  a  lattice  Green’s 
function.  The  method  is  fast  because  the  Green’s  function  is  sparse  in  reciprocal  space.  A  2D  fast 
Fourier  transform  (FFT)  of  atomic  displacements  in  the  plane  is  performed.  Forces  can  then  be 
calculated  by  multiplying  by  a  precalculated  Greens  function  for  each  wave  vector  q  and 
performing  an  FFT  back  to  real  space.  The  Greens  function  at  each  q  is  a  matrix  of  dimension  3n 
where  n  is  the  number  of  layers  that  interact  with  the  atomistic  region  (typically  1-3).  The  total 
time  scales  as  L"lnL  instead  of  L  . 

Our  implementation  of  the  method  differs  in  several  important  ways.  One  corrects  an  error 
that  we  discovered  in  the  algorithm.  When  interactions  extend  beyond  the  nearest  layer  of  atoms, 
previous  implementations  did  not  correctly  treat  the  change  in  interactions  due  to  the  change  in 
layer  spacing.  This  can  be  a  significant  effect  in  certain  cases.  We  also  developed  an  analytic 
method  for  calculating  the  Greens  function  at  low  temperatures.  The  previous  implementation 
used  separate  simulations  to  find  G.  The  final  improvement  was  to  implement  the  algorithm  on 
graphics  processor  units  (GPUs).  This  allowed  us  to  treat  L  as  large  as  8192  atoms  or  total 
volumes  of  550  trillion  atoms  on  a  single  GPU.  The  comprehensive  studies  described  below 
would  not  have  been  possible  without  this  speedup. 

A  paper  describing  the  modifications  in  the  algorithm  is  currently  in  preparation  (Pastewka 
and  Robbins  2012a)  and  an  excerpt  is  included  in  Appendix  I.  We  will  make  the  implementation 
available  as  an  open  source  add-on  to  LAMMPS.  The  group  maintaining  the  files  developed  by 
Kong  et  al.  (2009)  has  agreed  to  update  the  files  to  include  our  improvements  and  we  are  also 
talking  to  the  developers  of  LAMMPS  (lammps.sandia.gov)  about  incorporating  it  in  future 
releases. 

2.  MOLECULAR  DYNAMICS  SIMULATIONS  OF  SINGLE  ASPERITY  CONTACT 

One  thrust  has  been  molecular  dynamics  (MD)  studies  of  contacts  between  a  single  spherical 
asperity  and  a  substrate,  or  between  two  flat  surfaces.  In  addition  to  providing  understanding  of 
the  nature  of  these  contacts  at  the  atomic  scale,  these  simulations  provide  tests  of  continuum 
models  of  contact  that  can  be  used  to  study  multi-asperity  contacts  on  larger  length  scales. 

Three  specific  studies  are  described  below.  The  first  examines  the  profound  affect  that  a 
single  layer  of  adsorbed  atoms  can  have  on  contact  and  friction.  The  second  shows  that  thermal 
fluctuations  lead  to  large  pressures  that  suppress  the  number  of  contacting  atoms  by  orders  of 
magnitude.  Finally  the  effect  of  contact  and  tip  size  on  the  friction  coefficient  is  evaluated. 


Adsorbed  Monolayers 

One  study  (Cheng  et  al.  2010)  examined  the  effect  of  adsorbed  monolayers  on  a  spherical  tip 
with  a  radius  R  that  is  typical  of  atomic  force  microscope  (AFM)  experiments.  Previous  studies 
of  bare  surfaces  showed  that  the  atomic  scale  deviations  from  a  sphere  that  are  present  on  any  tip 
constructed  from  discrete  atoms  lead  to  significant  deviations  from  continuum  theory  and 
dramatic  variability  in  friction  forces  (Luan  and  Robbins  2005,  2006).  However  surfaces  are 
rarely  bare.  Most  are  coated  by  a  monolayer  of  water,  hydrocarbon  or  other  adsorbed  material 
after  any  exposure  to  air. 

We  found  that  introducing  an  adsorbed  monolayer  leads  to  larger  deviations  from  continuum 
theory,  but  decreases  the  variations  between  tips  with  different  atomic  structure.  Although  the 
film  is  fluid,  it  remains  in  the  contact  and  behaves  qualitatively  like  a  thin  elastic  coating  except 
for  certain  tips  at  high  loads.  Some  specific  findings  that  are  detailed  in  the  paper  are: 

1)  The  contact  area  increases  dramatically  with  a  single  molecular  layer  (Fig.  1).  The  area  at 
the  lowest  load  studied  was  higher  than  the  highest  area  predicted  by  continuum  theory. 

2)  The  contact  stiffness  is  reduced  by  an  order  of  magnitude  for  small  displacements,  but 
approaches  the  continuum  result  at  large  displacements.  The  extra  compliance  comes 
from  compression  of  the  monolayer  but  is  equivalent  to  continuum  results  for  a  thicker 
layer  that  is  much  softer  than  the  solid. 

3)  Monolayers  naturally  lead  to  a  friction  that  scales  linearly  with  load,  much  as  was  seen 
for  flat  surfaces  (He  et  al.  1999,  Muser  et  al.  2001).  In  contrast,  bare  surfaces  exhibit  very 
different  scaling  with  load  and  friction  forces  that  vary  by  two  orders  of  magnitude  (Luan 
and  Robbins  2005,  2006). 
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Figure  1.  Squares  show  contacting  atoms 
(a)  with  an  adsorbed  monolayer  and  (b) 
without  a  monolayer.  The  load  and  all 
other  parameters  are  the  same  and 
cr~0.3nm  is  an  atomic  diameter. 


Figure  2.  Displacement  S  of  tip  (closed  symbols)  and 
substrate  (open  symbols)  for  different  tip  geometries  and 
continuum  prediction  (solid  line).  The  much  larger  tip 
displacement  indicates  a  low  stiffness  from  compression  of 
the  adsorbed  layer.  R  is  the  tip  radius,  N  the  load  and  E  the 
effective  contact  modulus. 


Effect  of  Thermal  Fluctuations  on  Contact 

Thermal  fluctuations  are  generally  ignored  in  continuum  treatments  of  contact  mechanics 
and  in  many  simulations  that  attempt  to  test  continuum  theories  at  atomic  scales.  A  detailed 
study  (Cheng  and  Robbins  2010b)  of  the  effect  of  thermal  fluctuations  on  contact  was  completed 
for  the  simplest  possible  geometry  -  two  atomically  flat  surfaces.  The  results  show  that  it  is  very 
difficult  to  extend  the  continuum  definition  of  contact  area  to  atomic  scales.  The  results  were 
then  shown  to  describe  single  asperity  contact  of  a  sphere  on  a  plane. 

The  flat  surfaces  were  (001)  surfaces  of  fee  crystals  with  the  same  (commensurate)  or 
different  (incommensurate)  lattice  constants.  An  amorphous  surface  obtained  by  removing  all 
atoms  above  a  reference  plane  was  also  considered.  From  a  continuum  perspective  the  true  area 
of  contact  between  such  atomically-flat,  nonadhesive,  surfaces  should  be  equal  to  the  apparent 
area  at  any  positive  load.  However,  simulations  show  that  the  fraction  of  atoms  experiencing  a 
contact  force  at  any  instant  in  time  is  orders  of  magnitude  smaller  than  unity  and  rises  roughly 
linearly  with  load  over  several  decades  in  load. 

Fig.  3  illustrates  results  for  the  three  types  of  surfaces.  The  x-axis  corresponds  to  the 
compressive  strain  perpendicular  to  the  surface  and  extends  up  through  the  range  associated  with 
the  hardness  of  typical  materials.  Yield  usually  occurs  at  strains  of  10'4  to  1 0"2  and  the  fraction  of 
area  in  contact  rises  linearly  from  about  0.3%  to  30%  over  this  range.  Even  the  fraction  of  atoms 
that  contact  at  any  instant  of  the  entire  simulation  (triangles)  remains  well  below  unity  except  in 
the  very  special  case  of  two  identical  flat  crystals  (Fig.  3a). 

These  counterintuitive  results  are  quantitatively  explained  by  a  very  simple  model  for 
thennal  vibrations  (solid  lines  in  Fig.  3  and  4).  Each  atom  feels  a  potential  well  from  its  surface 
that  can  be  approximated  by  a  single  spring  keff.  From  dimensional  analysis,  this  must  be 
keff=cAaE  /a  where  a  is  the  molecular  diameter,  Aa  the  surface  area  per  atom,  E  the  contact 
modulus  and  c  a  constant  of  order  unity.  The  simulations  showed  c=0.92  for  Lennard-Jones 
interactions.  The  thennal  pressure  associated  with  the  rms  thermal  fluctuation  xrm,  about  the 
lattice  site  is: 

PT-keff  Xrms/Aa  —  CE*  (Xnns/cj). 

From  the  Lindemann  criterion,  xrms/<j  is  10%  at  the  melting  temperature  Tm,  and  rises  as  T  . 
Thus  pT/E*  ~  0. 1  (T/Tm)° 5  and  the  fraction  of  area  required  to  support  a  normal  load  L  is 

(i) 

This  simple  argument  is  in  reasonable  agreement  with  Fig.  3  where  T/Tm  ~  1/4. 

The  analytic  argument  produces  two  other  predictions  that  are  confinned  by  the  simulations. 
One  is  that  the  distribution  of  local  forces  at  any  instant  has  a  simple  exponential  form.  The 
second  is  that  the  time  averaged  force  on  an  atom  fc  is  directly  related  to  the  fraction  of  time  pc 
that  it  is  in  contact.  As  illustrated  in  Fig.  4,  this  relation  is  independent  of  whether  the  atom  is  on 
a  commensurate,  incommensurate  or  amorphous  surface  and  whether  the  surface  is  flat  or 
spherical.  Note  that  the  force  must  be  close  to  the  ideal  yield  stress  before  atoms  spend  more  than 
half  their  time  in  contact  and  the  actual  yield  stress  is  typically  orders  of  magnitude  smaller. 
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Figure  3.  Fraction  of  area  in  contact  Ac/A0  vs. 
mean  normal  strain  L/A0E*  for  (a)  commensurate, 
(b)  incommensurate  and  (c)  amorphous  surfaces. 
Circles  indicate  instantaneous  contact,  triangles 
show  contact  after  long  times  and  lines  are  a 
simple  analytic  theory.  L  is  the  load  and  E*  the 
contact  modulus 
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Figure  4.  Relation  between  time  averaged  force  fc  and 
fraction  of  time  in  contact  pc  for  flat  commensurate  (+), 
flat  incommensurate  (x),  flat  amorphous  (♦),  spherical 
commensurate  (A),  spherical  incommensurate  (  )  and 
amorphous  spherical  (o)  surfaces.  Here  T/Tm~l/4.  The 
solid  line  is  from  the  analytic  theory  developed  in  Cheng 
and  Robbins  (2010b).The  inset  shows  a  linear  plot  for  high 
pc.  The  ideal  yield  stress  corresponds  to  ~5  s/a,  where  e 
and  ct  are  the  binding  energy  and  atomic  diameter 


Scale  Dependence  of  Friction 

Traditional  continuum  theories  of  single  asperity  contact  have  no  intrinsic  microscopic  length 
scale.  All  results  depend  only  on  ratios  of  macroscopic  lengths  such  as  the  radius  of  curvature  R 
and  contact  radius  a  (Johnson  1985).  Hurtado  and  Kim  (1999a,  1999b)  suggested  that  friction 
should  be  scale  dependent  in  certain  limits.  They  argued  that  contacts  could  advance  through 
propagation  of  dislocations  across  the  interface  rather  than  unifonn  sliding.  The  Burgers  vector 
of  the  dislocation  (~c)  and  core  size  of  the  dislocation  introduce  microscopic  length  scales.  Just 
as  dislocations  reduce  the  yield  stress  in  bulk  crystals,  dislocations  led  to  a  reduced,  scale 
dependent  friction. 

In  addition  to  lowering  friction,  Hurtado  and  Kim’s  model  provided  a  scale  dependent 
mechanism  for  partial  slip,  where  sliding  occurs  within  the  outer  ring  of  a  contact  and  gradually 
spreads  inwards.  The  previous  Cattaneo-Mindlin  theory  for  partial  slip  does  not  have  a 
microscopic  length  scale  (Johnson  1985).  Improved  understanding  of  partial  slip  was  a 
motivation  of  the  proposed  research,  because  it  is  a  major  source  of  damping  in  mechanical 
joints,  including  mounts  for  aircraft  engines  (Segalman  2006). 

Testing  the  scale  dependence  proposed  by  Hurtado  and  Kim  poses  significant  challenges. 
They  predicted  that  the  friction  would  drop  roughly  as  (a/  g)'1/2  over  scales  from  a /  a  ~102  to  104. 
Since  R  must  be  bigger  than  10a  to  be  in  the  linear  elastic  limit,  this  means  linear  dimensions  of 
10  to  10  a  are  required.  Gao  (2010)  extended  Hurtado  and  Kim’s  work  and  found  transitions  at 
length  scales  about  an  order  of  magnitude  smaller,  but  still  taxing  for  standard  atomistic 
simulations. 


Using  the  Green’s  function  method  described  above,  Ligneres  and  Robbins  (2012)  examined 
the  friction  between  spherical  tips  of  radius  R/a  between  50  and  6000  (~2pm).  The  results  are 
very  different  than  envisioned  by  any  of  the  previous  continuum  theories  and  are  currently  being 
prepared  for  publication.  The  main  findings  are: 

1)  The  mechanism  of  sliding  is  never  dislocation  or  partial  slip  as  envisioned  by  the 
Cattaneo-Mindlin,  Hurtado-Kim  and  Gao  models.  The  Peierls  barrier  for  advance  of  the 
dislocation  is  too  small  to  pin  dislocations  that  are  formed  in  response  to  shear  stress. 

2)  The  area  and  load  dependence  of  friction  forces  depends  on  the  precise  geometry  of  the 
tips. 

3)  For  incommensurate  tip/substrate  pairs  (different  lattice  constants),  the  friction  is  small, 
nonuniversal  and  rises  less  than  linearly  with  area  and  load. 

4)  The  commensurate  case  of  equal  lattice  constants  is  the  only  one  that  allows  for  the 
surface  dislocations  in  Hurtado  and  Kim’s  treatment.  Once  the  static  friction  is  exceeded 
the  entire  tip  advances  by  a  distance  of  order  a  because  dislocations  cannot  be  pinned. 

5)  The  friction  coefficient  is  not  only  a  function  of  a/a,  as  in  the  results  of  Hurtado  and  Kim 
and  of  Gao.  There  is  an  explicit  dependence  on  R  that  is  completely  absent  in  their 
models. 

6)  Fig.  5  shows  that  the  static  friction  coefficients  for  commensurate  systems  of  all  sizes 
collapse  onto  a  universal  curve  as  a  function  of  a/(Ro)  .  Detailed  analysis  shows  that  the 
local  slope  of  the  surface  (a/R)  plays  a  role  in  the  force  at  which  the  tip  can  advance. 


a/(Rcj)1/2 


Figure  5.  Friction  coefficient  p  as  a  function  of  (a2/Ra)1/2  for  tips  with  radius  R/cj  equal  to  300  (♦),  600  (■),  1200 
(A),  2400  (x)  and  4800  (*).  The  coefficient  of  friction  is  just  less  than  0.7  for  flat  surfaces  with  the  same  spacing 
and  interactions. 


MULTIASPERITY  CONTACT 

Macroscopic  surfaces  are  typically  rough  on  a  wide  range  of  scales,  often  from  less  than  a 
micrometer  to  more  than  a  millimeter.  The  presence  of  roughness  on  a  wide  range  of  length 
scales  has  profound  effects  on  contact  and  friction  between  experimental  surfaces.  Under  a  broad 
range  of  conditions,  the  area  of  intimate  contact  between  nonadhesive  rough  surfaces  Ac  is  orders 
of  magnitude  smaller  than  the  apparent  surface  area  Ao  and  rises  linearly  with  applied  load  L. 
This  proportionality  between  area  and  load,  and  an  assumption  of  a  constant  yield  stress  in 
contacts,  provides  the  most  common  explanation  (Johnson,  1985)  for  Amontons'  laws  that 
friction  is  proportional  to  load  and  independent  of  A0.  Because  Ac  is  small,  the  interfacial  region 
is  very  compliant.  In  applications  from  jet  engine  mounts  to  micro-electromechanical  systems 
(MEMS),  the  interfacial  compliance  can  significantly  reduce  the  stiffness  of  joints  formed  by 
pressing  two  components  together  (Segahnan  et  al.  2006).  Electrical  and  thermal  conductances 
are  also  strongly  affected. 

Dimensional  analysis  implies  that  the  linear  relation  between  area  and  load  has  the  fonn: 

Ac  =  k  L  /E*  h’^  *  (2) 

where  k  is  a  constant  of  order  unity,  the  contact  modulus  E  characterizes  the  elastic  response  of 
the  substrate,  and  the  rms  surface  slope  h  ’  rms  is  a  dimensionless  quantity  characterizing 
roughness.  Numerical  tests  show  k  has  a  constant  value  of  about  2  for  most  surfaces  (Hyun  et  al. 
2004).  One  of  the  goals  of  the  grant  was  to  develop  equations  like  Eq.  (2)  that  describe  the 
stiffness  of  contacts  and  test  them  with  numerical  calculations.  A  second  was  to  include  the 
effect  of  adhesion  on  contact  area. 

The  following  sections  describe  the  three  major  studies  undertaken  during  the  grant.  The 
first  calculated  the  spatial  correlations  between  contacting  regions  in  order  to  test  analytic 
theories  of  contact  area.  The  results  show  that  conventional  models  give  qualitatively  wrong 
correlations  because  they  neglect  asperity  interactions.  The  second  develops  and  tests  a  very 
simple  model  that  describes  the  stiffness  of  mechanical  contacts  between  rough,  nonadhesive 
surfaces.  The  final  study  has  resulted  in  a  very  simple  and  robust  picture  for  the  contact  area 
between  rough,  adhesive  surfaces.  The  model  explains  why  adhesion  between  most  macroscopic 
surfaces  is  so  weak,  and  provides  a  simple  dimensionless  number  that  predicts  when  adhesion 
will  be  strong. 

Statistics  of  Asperity  Heights 

A  wide  variety  of  natural  and  man-made  surfaces  exhibit  self-affine  fractal  roughness.  Examples 
include  common  machined  surfaces,  fracture  surfaces  and  surfaces  formed  by  deposition  or 
erosion  (Mandelbrot  1982,  Greenwood  1984,  Meakin  1998).  The  hallmark  of  self-affine  surfaces 
is  that  the  change  in  height  h  over  a  lateral  length  1  scales  as  a  power  law:  h  cc  1H,  where  H  is 
called  the  Hurst  exponent.  For  typical  experimental  surfaces,  H  ranges  from  0.5  to  0.8.  A  random 
walk  gives  H=0.5. 

Calculations  for  the  contact  geometry  have  traditionally  been  based  on  the  bearing  area 
model  (Johnson  1985).  The  problem  of  contact  between  two  rough  surfaces  can  be  mapped  to 
contact  between  a  rough  and  flat  surface  (Johnson  1985).  At  a  given  separation  h0  between 
surfaces,  contact  is  assumed  to  occur  wherever  the  rough  surface  would  interpenetrate  the  flat 
surface.  This  corresponds  to  the  set  of  points  above  h0  on  the  rough  surface.  The  widely  used 
Greenwood-Williamson  (1966)  model  and  generalizations,  calculate  the  load  required  to  bring 
this  region  into  contact  without  including  the  affects  of  elastic  interactions  between  asperity 
peaks. 


Persson  (2001,2007)  has  recently  developed  a  scaling  theory  for  contact  that  includes 
asperity  interactions  in  an  approximate  manner.  These  interactions  can  be  important  because  the 
deformation  produced  by  a  force  on  an  asperity  decays  only  as  the  inverse  separation.  This 
means  that  nearby  areas  are  unlikely  to  make  the  next  contact. 

Despite  the  large  difference  between  the  Greenwood-Williamson  and  Persson  models  they 
give  similar  results  for  the  total  contact  area.  The  respective  values  of  k  are  about  2.5  and  1.6 
compared  to  the  numerical  result  of  2  (Hyun  et  al.  2004).  This  similarity  in  total  area  makes  it 
interesting  to  compare  the  spatial  correlations  in  the  contacting  regions  in  order  to  distinguish 
between  the  two  models. 

In  recent  papers  (Hyun  and  Robbins  2007,  Campana  et  al.  2008)  we  studied  the  square  of 
the  Fourier  transform  of  the  contacting  region,  which  gives  the  Fourier  transform  of  the  contact- 
contact  correlation  function  S(q).  Results  from  elastic  calculations  that  include  asperity 
interactions  yield  S(q)  oc  q'(  l+H)  for  a  wide  range  of  H.  Persson  later  showed  that  his  calculation 
gave  the  same  scaling  exponent  (Persson  2008). 

During  the  grant  we  calculated  the  correlation  function  S(q)  for  the  bearing  area  model  using 
numerical  results  and  an  analytic  theory  (Ramisetti  et  al.  2011).  Both  gave  S(q)  oc  q-(2+H>  raqlcr 
than  the  correct  result  for  interacting  asperities  S(q)  oc  q‘(1+H) .  This  failure  of  the  bearing  area 
model  has  important  physical  consequences.  For  example,  almost  all  area  in  the  bearing  area 
model  is  in  the  largest  connected  regions,  while  most  area  is  in  the  smallest  connected  regions  of 
the  full  elastic  calculation.  This  leads  to  very  different  results  for  measurable  quantities  such  as 
contact  conductance  of  heat  and  electricity,  as  well  as  very  different  predictions  for  adhesion. 

Mechanical  Stiffness  of  Multiasperity  Contacts 

The  component  of  the  stiffness  coming  from  the  interface  between  rough  solids  was 
calculated  for  surfaces  with  a  wide  range  of  Hurst  exponents  using  fully  atomistic  calculations 
and  the  Green’s  function  method.  A  simple  relation  like  Eq.  (2)  described  all  the  results  quite 
accurately  (Akarapu  et  al.  2011)  and  can  be  used  to  predict  the  response  of  experimental 
surfaces. 

The  interfacial  stiffness  k  is  found  to  scale  linearly  with  the  applied  load  L.  Since  k  is  the 
opposite  of  the  derivative  of  load  with  separation  u,  one  must  have  a  relation  like: 

k  =  -dL/du  =L/yhrms  (3) 

where  y  is  a  constant  of  order  unity  and  the  rms  height  hrms  is  a  natural  length  scale 
characterizing  the  surface.  This  type  of  relation  had  been  suggested  by  Berthoud  and 
Baumberger  (1997)  and  Persson  (2007).  Our  simulations  provide  the  first  detailed  numerical  test 
of  Eq.  (3).  As  shown  in  Fig.  6  this  relation  is  quite  general,  with  y  near  0.5.  Several  other  aspects 
of  the  results  are  more  surprising: 

1)  The  simulations  included  both  perfectly  elastic  systems  treated  with  Green’s  functions 
and  atomistic  systems  that  defonned  plastically,  yet  the  two  cases  gave  similar  values  of 
k  and  y.  A  detailed  analysis  indicates  that  plasticity  is  confined  to  regions  near  the 
surface.  It  is  important  in  eliminating  extremely  large  local  pressures  but  does  not  affect 
the  long-range  scaling  (Sharp  et  al.  2012). 

2)  Atomic  steps  affect  plasticity  but  have  much  less  influence  on  the  total  area  than  in  two 
dimensional  studies  (Luan  2009). 

3)  The  lateral  stiffness  k  [  associated  with  elastic  deformation  of  the  contacts  can  be  larger 
than  the  nonnal  stiffness  kN.  Continuum  theory  predicts  the  opposite,  but  does  not 
include  lattice  anisotropy. 


4)  The  stiffness  only  depends  on  the  lateral  position  of  contacting  regions  and  not  on  the 
height  variations  or  other  factors.  This  was  used  to  derive  a  linear  relation  between  area 
and  a  normalized  stiffness  k*  (Fig.  7). 

5)  There  is  an  additional  compliance  associated  with  relative  motion  of  contacting  atoms  at 
the  interface.  This  atomistic  compliance  has  almost  no  affect  on  the  normal  stiffness  (Fig. 
7a),  but  can  reduce  the  lateral  stiffness  and  friction  by  orders  of  magnitude  (Fig.  7b).  The 
results  are  analogous  to  single  asperity  results  of  Luan  and  Robbins  (2005,2006).  A 
scaling  study  indicates  that  atomistic  effects  can  be  an  important  factor  in  macroscopic 
samples  (Akarapu  et  al.  2011). 


Figure  6.  The  scaled  (a)  normal  stiffness  and  (b) 
transverse  stiffness  as  a  function  of  scaled  load.  Results 
are  for  H=0.5  (open  symbols)  and  H=0.8  (filled  symbols) 
in  systems  with  linear  dimension  of  189.2cr  (°),  94.6a  (□), 
or  47.3a  (A).  Dashed  lines  have  slope  1/y  with  y=0.48. 
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Figure  7.  The  scaled  (a)  normal  stiffness  and  (b) 
transverse  stiffness  as  a  function  of  Ac/A0.  Results  are 
for  H=0.5  (open  symbols)  and  H=0.8  (filled  symbols)  in 
systems  with  linear  dimension  of  189.2a  (°),  94.6a  (□), 
or  47.3a  (A).  Pluses  show  Green’s  function  results.  The 
total  tranverse  stiffness,  multiplied  by  20  to  make  it 
visible,  is  shown  for  189.2a  systems  with  H=0.5  (x)  or 
0.8  (hexagons)  and  a  triangular  lattice  with  H=0.8  (*). 
Dashed  lines  have  slope  0.87. 


Adhesive  Contact  Between  Rough  Surfaces 

Adhesive  van  der  Waals  interactions  are  present  between  almost  all  surfaces.  The  forces  required 
to  break  these  interactions  are  in  the  MPa  range,  and  enable  Geckos  to  climb  walls  (Autumn 
2006)  and  tape  or  putty  to  stick  to  surfaces.  Similar  forces  are  present  between  nearly  all  surfaces 
around  us,  and  yet  we  rarely  observe  adhesive  forces  when  we  lift  objects  off  of  tables.  Fuller 
and  Tabor  resolved  this  adhesion  paradox  by  showing  that  surface  roughness  prevents  adhesion 
by  limiting  the  area  that  can  come  into  contact  (Fuller  and  Tabor  1975).  However  a  detailed 
theory  for  this  suppression  and  a  criterion  for  the  presence  of  adhesion  for  soft  materials  has 
remained  lacking. 

We  have  used  the  Green’s  function  technique  described  above  to  explore  a  full  range  of 
parameter  space  with  different  rms  roughness  ho,  rms  slope  h’o,  Hurst  exponent  H,  adhesion 
energy  y,  adhesion  range  Ac  and  system  size.  The  results  guided  the  development  of  a  simple 


scaling  theory  that  captures  all  our  simulations  and  provides  a  simple  prediction  for  when 
adhesion  will  occur  and  how  the  contact  area  scales  with  load.  The  work  is  being  prepared  for 
publication  (Pastewka  and  Robbins,  2012b)  and  a  draft  is  provided  in  Appendix  II. 

One  important  finding  is  that  the  portion  of  the  surface  that  experiences  repulsive  forces, 
Arep,  is  always  fractal  and  has  a  characteristic  width  w  that  is  comparable  to  the  smallest  scale  of 
surface  roughness.  As  a  result,  the  perimeter  of  the  repulsive  region  P=pArep/w,  where  p  is  about 
2.8  for  all  surfaces  studied.  We  also  find  that  the  ratio  of  the  Arep  to  the  repulsive  load  on  the 
same  region  is  given  by  the  same  k  as  for  nonadhesive  surfaces. 

Given  the  repulsive  area,  continuum  theory  for  nonadhesive  surfaces  (Johnson  1985)  can  be 
used  to  estimate  the  range  R  away  from  the  perimeter  where  attractive  forces  will  be  felt.  The 
ratio  of  repulsive  to  attractive  area  Aatt=PR  should  then  be  given  by: 

Arey  _  2  //lpw\2^3 
Aatt  P  '  2AC  ) 


Figure  8  shows  that  results  for  all  surfaces  fall  very  close  to  this  prediction  with  P=2.8. 

The  total  attractive  contribution  to  the  load  is  Natt=AattY/Ac,  since  y/Ac  is  the  average 
adhesive  force  over  the  range  of  interaction.  One  can  define  a  renonnalized  k’  that  relates  the 
repulsive  area  to  the  total  load  N=Nrep-Natt 

,  N 

Arep  ~  K  h^E~* 

Then  plugging  in  the  relations  above  one  finds 

1  _  l  +  J_  _  l  P  8  /2AC\2/3 

k'  k  Katt  k  2h'Q  Ac  \h'0w)  [  : 

where  8  =  y/E  is  a  length  characterizing  the  relative  strength  of  adhesion  and  elasticity.  Fig.  9 
shows  a  test  of  the  predictions  of  this  simple  relation  for  a  wide  range  of  surfaces. 

Equation  5  allows  the  attractive  area  to  be  predicted  in  terms  of  a  few  simple  quantities. 
Adhesion  will  have  negligible  effect  when  the  term  on  the  far  right  is  much  smaller  than  1/k~0.5. 
When  the  right  hand  side  becomes  negative,  there  will  be  substantial  adhesion.  Simulations  in 
this  limit  show  a  finite  repulsive  area  at  zero  load,  and  a  negative  load  is  needed  to  separate  the 
surfaces. 

For  most  solids  y  and  E*  are  determined  by  the  same  interactions  and  y/E*  <  0.01a.  The 
range  of  interactions  is  typically  ~0.1  a  or  larger  and  the  width  w  >  a.  As  a  result,  the  correction 
from  adhesion  is  negligible  unless  the  rms  slope  is  extremely  small.  This  type  of  situation  may 
arise  in  wafer  bonding  of  silica  with  sub  nanometer  roughness,  but  not  for  most  of  the  surfaces 
around  us. 

Tape,  putty  and  elastomeric  materials  typically  have  anomalously  small  moduli  because  they 
contain  chain  molecules  that  are  relatively  fluid  in  between  a  dilute  set  of  chemical  crosslinks. 
The  modulus  reflects  the  entropic  cost  of  stretching  chains  between  crosslinks  while  y  is 
determined  by  van  der  Waals  interactions.  As  a  result,  y/E*  can  be  four  orders  of  magnitude 
larger  than  for  other  materials  and  Eq.  4  becomes  negative  even  for  typical  rough  surfaces. 


Figure  9.  Test  of  Equation  5  for  the  same  set  of  systems.  Note  that  as  kalt  approaches  -0.5,  and  k’ 
approaches  zero  there  is  a  systematic  deviation  towards  greater  adhesion. 
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APPENDICES 


I.  METHODOLOGY  FOR  CONTACT  MECHANICS  OF  CRYSTAL  LATTICES 

Analysis  of  the  response  of  an  isotropic  elastic  half-space  to  spherical  contact  geometries 
and  point  forces  goes  back  to  the  work  of  Hertz  [1] ,  Boussineq  and  Cerruti  in  the  19th  century. 
Since  then,  continuum  expressions  have  been  derived  for  various  contact  geometries  [2] ,  and 
have  become  indispensible  tools  in  mechanical  engineering.  Continuum  mechanical  contact 
theories,  however,  break  down  at  the  atomic  scale  [3].  Yet,  an  understanding  of  long-ranged 
clastic  deformations  of  the  bulk  down  to  the  atomic  scale  is  important,  for  example,  for  the 
interpretation  of  scanning  probe  experiments  or  the  contact  stiffness  of  atomically  rough 
surfaces  [4] .  Also,  cherni-  and  physisorption  processes  at  crystal  surfaces  can  be  appreciably 
influenced  by  the  elastic  response  of  the  supporting  solid. 

Here,  we  develop  a  transfer  matrix  formulation  of  the  Green’s  function  that  can  be  used  to 
renormalize  the  clastic  interaction  of  any  crystalline  solid  for  which  the  dynamical  matrix  is 
known.  Such  transfer  matrix  formulations  have  been  traditionally  used  in  the  analysis  of  the 
electronic  [5]  and  phononic  [6]  structure  of  surfaces,  but  have  to  the  best  of  our  knowledge 
not  been  discussed  in  the  context  of  interactions  that  reach  farther  than  nearest  neighbors. 
Such  long-ranged  interactions  are  necessary  for  a  realistic  description  of  materials.  Hence, 
the  treatment  of  long-ranged  interactions  allows  to  employ  the  renormalized  half  space  as 
an  ideally  matched  clastic  boundary  in  atomistic  calculations. 

A.  Elastic  surface  Green’s  functions 

1.  Green’s  function  of  the  full  dynamical  system 

Let  us  expand  vibrations  of  the  crystalline  solid  around  the  ground-state  configuration 
r^a  of  each  atom  ia.  To  obtain  the  usual  harmonic  approximation,  we  write  all  dynamical 
equations  in  terms  of  the  the  displacement  uia  of  atom  ia  from  its  ground-state  position. 
Here  and  in  the  following,  atoms  are  indexed  by  ia  where  a  is  a  layer  index  (as  counted 
from  the  surface  of  the  crystal  where  a  =  0)  and  i  denotes  inequivalent  atoms  within  each 
layer.  The  full  dynamical  Green’s  function  can  then  be  obtained  from  the  Fourier  transform 
of  the  dynamical  equation 

ViojpUjp  =  fia ■  (1) 

iP 
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Here,  m  is  the  mass  of  each  atom  and  we  limit  ourselves  to  the  treatment  of  a  single 
atomic  species.  A  generalization  of  the  following  procedure  to  multi-species  unit  cells  is 
straight  forward . 

We  now  treat  the  dynamical  equation  in  reciprocal  space  parallel  to  the  surface  and  in 
real  space  perpendicular.  With 


U io  (t) 


d3q  du 
(27 r)3  27 r 


u(q)e 


iq-rioc+icut 


(2) 


where  q  lies  in  the  plance  of  the  surface  and  flo  =  +  UiQ  is  the  position  of  atom  ia,  we 
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with  Fia  =  J2  D ia  jpUjp  and  Fjpia  =  Fjp  —  ria .  From  Eq.  (4)  we  identify 
30 


D  ap{q)  =  Y,  D  ianpe1^*  (5) 

n^Eneighb. 

where  n  run  over  those  neighbors  in  layer  (3  where  D ianp  ^  0  and  rianp  is  the  distance  vector 
to  those  neighbors.  In  combination  with  the  dynamical  Eq.  (1)  this  yields  the  eigenvalue 
equation 

(-mu2  +  D (q))  u(q,  u)  =  f(q,  u)  (6) 

for  the  permissible  modes,  and  the  equation 


1  +  mu2G(q.u)  =  D(g)G(<f,  u;) 


(7) 


for  the  dynamical  Green’s  function  G(q,  u). 

In  the  remainder  of  this  paper  we  concern  ourselves  with  the  static  response  only.  The 
static  response  is  given  by  the  limit  u  =  0,  such  that  the  governing  equation  becomes 


G(q)G(q)  =  1, 


(8) 


where  we  dropped  the  frequency  variable  u  for  brevity.  In  order  to  understand  the  displace¬ 
ments  u0(q)  of  the  surface  layer  due  to  an  external  force  fo(d)  we  only  need  the  Green’s 
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function  Goo  (q)  element  that  couples  the  surface  layer  to  itself.  It  is  usually  beneficial  to 
treat  the  surface  stiffness  $>(q)  =  o(q),  such  that 

fo(q)  =  &(g)uo(?)>  (9) 


instead.  Eq.  (9)  resembles  Hooke’s  law,  and  the  coefficients  $  can  be  regarded  as  renor¬ 
malized  spring  constants  that  govern  the  response  of  the  elastic  half  space.  Note  that  even 
though  the  atomistic  interaction  within  the  bulk  may  be  short  ranged,  real  space  $  coeffi¬ 
cients  typically  couple  the  surface  over  all  length  scales. 

In  the  following  we  assume  that  the  total  crystal  consists  of  N  +  1  identical  crystal 
layers.  Then  U(g)  =  Diana(gJ,  the  dynamical  matrix  that  couples  within  each  layer,  V (q)  — 
Djan(a+i),  the  matrix  coupling  to  the  nearest  layers  beneath,  and  Vt(g),  the  matrix  coupling 
to  the  nearest  layer  above,  do  not  depend  on  the  index  a  for  a  >  0.  For  a  =  0  we  introduce 
Uo(/7)  =  Djono  and  note  that  V0  =  V.  Here  f  denotes  the  Hcrmitian  conjugate. 

The  dynamical  matrix  then  assumes  the  form 
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(10) 


where  we  dropped  the  explicit  q  dependence  for  brevity.  For  simplicity,  we  assume  no 
coupling  between  next-nearest  layers,  and  later  show  that  next-nearest  layer  interactions 
can  be  straightforwadly  mapped  on  a  theory  that  only  treats  nearest  layers.  We  now  give 
two  equivalent  formulations  that  lead  to  the  renormalized  surface  Green’s  function. 
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2.  Displacement  transfer  matrix  formulation 


The  displacement  transfer  formulation  starts  from  Eq.  (6)  with  uj  =  0.  We  can  write  this 
equation  explicitly  as 

U0m0  +  Vmi  =  /  (11) 

VW0  +  Uiq  +  Vm2  =  0 
V*«l  +  U-U2  +  V«3  =  0 


un— i  +  U un  T  V u^+i  0 


(12) 


V^Ujv-i  +  U  nun  —  0 

It  is  straightforward  to  see  that  this  set  of  equations  can  be  expressed  as 
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(13) 


(14) 


(15) 


where  T  is  called  the  transfer  matrix.  We  impose  the  boundary  condition  Mv  =  0,  yielding 

Mi  =  T12MV-1  and  m0  =  T22UN- 1-  (16) 


Inserting  into  Eq.  (11)  allows  to  eliminate  Mat_i  giving 

^(U0T22  +  VTi2)T2-21  (17) 

the  effective  response  of  the  surface  layer  to  an  external  force  /  under  conditions  of  a  fixed 
bottom  boundary. 
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3.  Green’s  function  transfer  matrix  formulation 


The  Green’s  function  formulation  starts  from  Eq.  (8).  Unlike  the  dynamical  matrix,  the 
Green’s  function  is  not  sparse.  We  denote  the  individual  elements  as 


0 

o 

o 

Goi 

Go2  •  • 

0 

o 

G10 

Gn 

g12  •• 

■  GiAT 

g20 

g21 

g22  •• 

■  G2V 

Gvo 

Gao 

Gat2  '  ' 

■  Gnn 

and  are  now  interested  in  extracting  the  surface  Green’s  function,  Goo-  Note  that  Gnn  and 
Goat  can  be  determined  equivalently. 

From  Eq.  (8)  we  obtain  generally  (N  +  l)2  equations  for  our  finite  system  with  N  +  1 
layers.  We  will  only  pick  the  N  +  1  equations  involving  the  surface  layer.  Those  are 

UoGoo  +  VGio  =  l  (19) 

VtG0o  +  UG10  +  VG2o  =  0 
VtG10  +  UG2o  +  VG3o  =  0 


VtGn_li0  +  UGtt,0  +  VGn+1)0  =  0  (20) 

V^Gat-1,0  +  UatGat,o  =  0  (21) 

We  now  define  the  transfer  matrix  Tn  as 

Gn+i,o  =  TnGn>o-  (22) 

The  surface  Green’s  function  Goo  and  stiffness  $  is  then  obtained  from  Eq.  (19)  as 

$  =  G001  =  U0  +  VT0.  (23) 

Combining  Eqs.  (20)  and  (22)  yields 

VT„Tn_1  +  UT„_1+Vt  =  0  (24) 

from  which  we  obtain 

VT„_!  =  -V(U  +  VTn)_1Vt.  (25) 
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This  continued  fraction  can  be  terminated  by  Eq.  (21)  which  gives 

VTV_!  =  -VU^Vt.  (26) 

Hence  we  can  iterate  Eq.  (25)  starting  from  Eq.  (26)  to  obtain  $  from  Eq.  (23). 


4 ■  Renormalized  dynamical  matrices  for  long-ranged  interactions 


If  the  interatomic  interaction  acts  over  a  distance  of  two  layers,  the  dynamical  matrix 


becomes 
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(27) 


Hence  we  can  renormalize  the  dynamical  matrices  U  and  V  according  to 


U0  = 


,U  = 


and  V  = 


(28) 


to  yield  an  effective  nearest  layer  model  with  more  that  three  degrees  of  freedom  per  unit 
cell. 

For  interactions  that  extend  over  three  neighboring  layers  we  get 
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which  leads  to  the  renormalization 


U0  = 

<  d 

— o 

d  < 

<  3 

,u  = 

1  U  V 

vt  u  v 

and  V  = 

0  0^ 
W  X  0 

(30) 

\ w t  vt  u 2J 

yWt  Vt  U  j 

vV  W  X j 

For  interactions  extending  over  even  farther  distances  the  procedure  can  be  carried  out 
accordingly. 

5.  Remarks 

The  displacement  and  Green’s  function  transfer  matrix  formalisms  of  sections  IA2  and 
I A  3,  respectively,  give  identical  results.  However,  at  the  Brillouin  zone  boundary  the  motion 
of  individual  atoms  decouples  as  the  Einstein  oscillator  limit  is  reached.  Hence  we  have 
V  =  0  and  the  transfer  matrix  T  as  defined  in  Eq.  (14)  becomes  singular.  For  practical 
purposes  the  Green’s  function  formulation  of  section  I  A3  is  hence  better  suited  and  will  be 
used  in  what  follows. 

Both  algorithm  are  iterated  exactly  N  times.  Hence,  the  perpendicular  dimension  of  the 
elastic  half-space,  i.e.  the  number  of  layers  N ,  is  explicitly  included  in  the  formulation. 
However,  the  quantity  <f>  is  well  defined  in  the  limit  of  an  infinite  elastic  half-space  N  — >  oo. 
Slowest  convergence  is  obtained  for  the  T  point  which  converges  to  <3>(r)  — >  0  as  IV  — »  oo. 
In  this  case,  the  half-space  becomes  infinitely  compliant  with  respect  to  a  uniform  force. 

B.  Stiffness  of  the  SC  (100)  and  FCC  (100)  surfaces 

We  now  use  the  Green’s  function  transfer  matrix  formulation  to  compute  the  stiffness 
matrix  for  the  simple  cubic  (SC)  (100)  and  the  face-centered  cube  (FCC)  (100)  surfaces. 
The  SC  lattice  is  stabilized  by  nearest  and  next-nearest  neighbor  interactions,  both  with 
spring  constants  k  —  1.  This  lattice  is  isotropic,  and  an  analytical  expression  that  agrees 
with  our  numerical  results  has  been  derived  in  Ref.  7.  The  FCC  lattice  interacts  via  nearest 
neighbors  with  spring  constant  k  —  1  only. 

For  completeness,  we  compare  the  stiffness  coefficients  to  the  behavior  of  <fr(q)  as  expected 
for  a  continuum  theory.  The  solution  is  obtained  by  a  Fourier-transform  of  the  Boussinesq 
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q  q  q 


FIG.  1:  Stiffness  matrix  coefficients  for  the  FCC  (100)  (blue),  the  SC  (100)  (green)  surface  and  an 
isotropic  elastic  solid  (red). 


solution  for  a  concentrated  point  force  on  an  isotropic  elastic  half-space  [2],  and  reads 
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2  q2 

q 

\ 
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(31) 


where  \i  is  the  shear  modulus  and  v  Poisson’s  ratio.  Figure  1  shows  the  value  of  the  diagonal 
and  off-diagonal  components  of  the  $  matrix  for  the  three  models. 


C.  Elastic  boundary  conditions  for  atomistic  calculations 

While  the  linear  elastic  response  of  an  atomic  half-space  is  important  to  study  contact 
mechanical  problems  in  the  linear  response  regime,  such  an  approach  can  easily  be  used 
as  an  elastostatic  boundary  condition  in  molecular  statics  or  dynamics  simulations  of  for 
examples,  indentation,  crack  propagation  or  nonlinear  contact  mechanics.  The  advantage 
of  the  scheme  presented  here  over  the  fluctuation  formulation  of  Campana  and  Muser  [8]  is 
that  the  Green’s  function  is  obtained  exactly  for  the  interatomic  potential  under  consider¬ 
ation  hence  minimizing  mismatch  of  the  elastic  response  at  the  interface.  Furthermore,  the 
response  close  to  the  T  point  is  only  poorly  sampled  by  an  evaluation  using  the  fluctuation 
theorem,  but  becomes  important  in  particular  if  the  alteration  in  the  atomistic  part  of  the 
system  reaches  out  beyond  nearest  neighbors  and  we  need  to  compute  the  Green’s  function 
for  multi-atom  unit  cells  as  discussed  in  section  I A  4. 

In  order  to  properly  couple  an  atomistic  system  with  interactions  that  reach  farther  than 
nearest  neighbors  to  the  elastic  half-space  we  need  one  final  consideration.  The  surface 
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FIG.  2:  For  balance  for  renormalized  surface  with  second-nearest  neighbor  interactions  and  explicit 
atoms  on  top. 

stiffness  <f>  is  based  on  the  harmonic  approximation  in  that  the  total  energy  of  the  system 
is  expanded  up  to  second  order.  In  the  bulk,  first  order  terms  cancel  since  the  system  in 
equilibrium  (see  also  Fig.  2).  This  does  not  hold  for  the  surface,  where  terms  linear  in  the 
displacements  give  rise  to  surface  relaxation  phenomena.  Note  that  such  relaxation  cannot 
occur  for  crystals  that  interact  via  nearest-neighbors  only,  since  in  this  case  the  neareast 
neighbor  distance  of  the  crystal  equilibrium  coincides  with  the  minimum  of  the  interaction 
potential  and  hence  forces  are  zero.  In  particular  when  adding  a  layer  of  explicit  atoms  on 
top,  the  force  due  to  these  atoms  needs  to  be  cancelled  on  each  renormalized  atom. 

The  situation  is  depicted  in  Fig.  2.  The  green  arrow  shows  the  force  that  the  explicit 
atom  excerts  on  one  of  the  renormalized  atoms.  In  the  bulk,  this  force  would  be  cancelled  by 
the  red  arrows,  which  are  however  missing  in  the  elastic  manifold.  We  hence  need  to  enforce 
cancellation  by  adding  a  constant  force  of  the  same  magnitude  and  opposite  direction  to  the 
renormalized  atom. 

II.  CONTACT  MECHANICS  OF  ADHESIVE,  RANDOMLY  ROUGH  SURFACES 

The  weak  van-der-Waals  interaction  allows  the  Gecko  to  overcome  gravity  and  stick  to 
walls  and  ceilings  [9-12],  The  same  force  acts  between  everyday  objects,  and  yet  we  are 
able  to  lift  pens,  plates,  and  almost  anything  else  with  ease.  Fuller  and  Tabor  resolved  this 
adhesion  paradox  when  finding  that  the  apparent  lack  of  macroscopic  adhesion  (or  stickiness) 
is  due  to  most  surfaces’  intrinsic  roughness  [13].  Yet,  a  detailed  microscopic  understanding  of 
how  the  non-adhesive  limit  is  controlled  by  the  surface’s  geometry  is  lacking.  In  particular, 
it  is  unknown  when  a  cross-over  to  sticky  behavior  happens,  and  how  this  cross-over  depends 
on  the  statistical  properties  of  the  surface  and  the  microscopic  interaction.  Here,  we  derive  a 
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simple  scaling  relation  that  describes  the  crossover  from  macroscopic  non-adhesive  to  sticky 
behavior.  This  crossover  is  controlled  by  the  characteristics  of  the  rough  surface  on  the 
smallest  scale  only.  Such  insights  are  important  for  the  design  of  adhesives  [9],  but  also 
to  understand  and  prevent  failure  at  interfaces,  such  as  debonding  of  thin  films  and  of 
composite  materials  [14],  or  fretting  and  pitting  in  tribological  systems. 

The  microscopic  origins  of  adhesion  are  manifold  and  range  from  simple  intermolecular 
interactions  [15]  over  capillary  bridges  [16]  and  adhesive  films  [17]  to  the  influence  of  bubbles 
of  air  at  the  contact  [18].  Microscale  adhesion  is  characterized  by  the  energy  per  surface  area 
7  that  is  gained  when  two  nominally  flat  surfaces  are  brought  into  contact.  Such  interactions 
give  rise  to  macroscale  adhesion,  or  stickiness,  that  manifests  itself  in  nonzero  pulloff  forces 
for  the  contacting  surfaces. 

Cohesive  or  Dugdale  zone  models  are  used  to  capture  the  effect  of  the  adhesive  energy  7 
in  continuum  theories  [19].  An  important  example  is  the  Maugis-Dugdale  (MD)  model  of  the 
adhesion  of  a  sphere  on  an  elastic  substrate  [20]  that  has  been  extensively  used  to  interpret 
scanning  probe  adhesion  experiments  [21].  Many  real  surfaces  are  not  flat  as  assumed  in 
continuum  theories.  Roughness  at  the  atomic  length  scale  can  have  appreciable  influence  on 
contact  and  adhesion  of  objects  as  simple  as  spherical  tips  with  atomic  scale  diameters  [3,  22], 
The  roughness  on  macroscopic  surfaces  typically  extends  over  many  length  scales  and  is 
suitably  approximated  by  Gaussian  self-affine  fractals  where  the  root  mean  square  height 
fluctuation  Sh  between  two  points  at  distance  A  on  the  surface  scales  as  5h  ~  \H  for  A 
between  As  and  A l  [23].  We  here  investigate  the  adhesion  of  such  surfaces  that  are  held  rigid 
on  a  square  flat  deformable  ideally  clastic  substrate  for  wide  variations  in  Hurst  exponent 
H,  root  mean  square  slope  h'0,  small  wavelength  cutoffs  \s  ,  and  large  wavelength  cutoff  A l- 
For  the  interaction  between  rigid  surface  and  clastic  substrate  we  choose  a  Lennard-Jones 
type  potential  as  an  approximation  for  the  van-der-Waals  interaction  (see  methods).  The 
sites  on  the  clastic  substrate  are  arranged  on  a  square  grid  with  lattice  constant  a 0  and  will 
be  regarded  as  atoms. 

The  clastic  interaction  is  solved  using  a  Fourier-transfornr  technique  [8]  with  a  linearized 
surface  Green’s  function.  We  only  treat  displacements  in  ^-direction  normal  to  the  surface 
such  that  the  total  elastic  energy  is  given  by  Ee]  =  A*  ^2q<lE*ul(q)  where  A0  is  the  total 
area  of  the  elastic  half  space  and  q  is  a  reciprocal  space  vector  in  the  plane  of  the  half  space. 
This  corresponds  to  the  isotropic  elastic  stiffness  matrix  Eq.  (31)  with  only  a  single  degree 
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of  freedom.  Each  discrete  point  on  the  clastic  manifold  will  be  regarded  as  an  “atom” . 
Each  atom  occupies  unit  area,  and  interacts  with  the  rigid  surface  through  the  potential 
V(x,y,z)  =  v(z  —  h(x,y ))  where  h(x,y )  is  the  rough  surface.  For  the  interaction  of  two 
flat  surfaces,  the  adhesive  energy  w  per  unit  area  is  given  by  the  minimum  of  the  atomic 
interaction  v(z).  Rigid  Gaussian  random  surfaces  with  Hurst  exponents  H  and  rms  slopes 
K  =  AFW)  were  generated  using  a  Fourier-hltering  algorithm.  All  wavevectors  with 
q  >  2-k/Xs  have  zero  power. 

After  hireling  the  equilibrium  contact  geometries  for  different  external  loads,  we  partition 
the  sites  on  the  elastic  surface  into  regions  of  atoms  that  feel  a  repulsive  and  regions  that  feel 
an  attractive  force,  as  well  as  regions  that  do  not  interact.  Fig.  3a  shows  these  individual 
regions  at  a  fractional  repulsive  area  of  5%.  Let  us  first  characterize  the  geometry  of  the 
repulsive  patches  by  determining  the  length  P  of  their  perimeter.  First,  we  note  that 
Prep  ~  Arep  in  stark  contrast  to  the  slit  island  analysis  which  yields  P  ~  A~H  [23].  The 
latter  is  a  consequence  of  the  fact  that  the  area  of  the  island  has  dimension  2  while  the 
perimeter  (or  “coastline”)  has  dimension  2  —  H.  This  implies  that  for  the  clastic  contact 
problem  the  total  repulsive  patch  is  a  fractal  object.  To  find  the  length  scale  that  controls 
A/P  consider  the  distribution  p(w)  of  the  lengths  w  of  segments  that  are  obtained  from  one¬ 
dimensional  slices  through  the  repulsive  contact  geometry  (see  Fig.  3b).  The  total  repulsive 
area  is  then  given  by  Arep  —  a0  f  dwwp(w)  =  a0w,  and  each  segment  contributes  2o0  to 
the  perimeter,  hence  Prep/Arep  =  ol/w  with  a  —  2.  Numerically,  we  find  a  slightly  different 
value  of  a  ~  2.8  which  is  attributable  to  segments  of  the  perimeter  that  area  parallel  to  the 
w-segments  missing  from  the  total  estimate.  This  is  an  effect  of  the  discrete  character  of 
the  contact  patches. 

Remarkably,  we  find  that  the  length  w  only  weakly  depends  on  the  contact  area  and  ad¬ 
hesion  strength  (see  Fig.  3c).  This  observation  holds  for  all  surfaces  studied  here.  Although 
the  contacting  surface  has  roughness  on  all  lengths,  only  a  single  length  scale  w  that  is  of 
the  order  of  the  short  cutoff  appears  to  be  important  for  a  description  of  the  contact 
problem. 

Besides  giving  the  ratio  between  area  and  circumference,  the  length  w  plays  the  role  of 
a  contact  diameter.  From  continuum  considerations,  it  is  straightforward  to  show  that  the 

3 

air  gap  A  between  the  two  separating  surfaces  raises  asymptotically  as  A  =  g(R)  ~  R 5 
for  small  distances  R  from  the  contact  edge.  Numerically,  we  find  exactly  this  power-law 


12 


APPENDICES 


FIG.  3:  (a)  Map  of  the  contacting  atoms  that  feel  a  repulsive  (orange)  and  attractive  (black)  force 
when  a  rough  surface  is  pushed  against  a  flat  elastic  manifold  (the  surface  shown  has  H  =  0.8, 
H'q  =  0.1,  As  =  4  and  A l  =  2048).  Atoms  that  feel  a  repulsive  force  cover  an  area  of  roughly  5%. 
(b)  Zoom  into  a  64  x  64  region  of  that  surface.  The  blue  vertical  lines  show  individual  segments  of 
length  w  (see  text),  (c)  Dependence  of  the  average  segment  length  w  on  contact  area,  (d)  Shape  of 
the  gap  between  rigid  and  elastic  surfaces  g  as  we  move  away  from  the  edge  of  the  contact  patch  in 
the  non-adhesive  hard-wall  limit.  The  inset  shows  the  behavior  for  R  <  0.05u)  close  to  the  contact 
edge.  The  dashed  line  is  the  gap  obtained  for  a  2D  wedge  indenter,  and  the  dashed-dotted  line 
shows  a  3/2  power  law. 

behavior  for  the  rough  contact  for  small  distances  from  the  contact  edge  with  a  crossover 
to  linear  behavior  at  larger  distances.  Indeed,  the  gap  data  collapses  as  we  normalize  all 
distances  by  the  contact  radius  w/2  and  all  heights  by  h'0w/2  as  shown  in  Fig.  3d. 

The  partitioning  into  repulsive  and  attractive  regions  of  the  area  A  =  Arep  +  Aatt  entails 
an  indentical  partitioning  of  the  total  load  N  =  Nrep  —  Natt.  We  find  that  the  repulsive  load 
Nrep  shows  a  linear  Amontons-like  dependence  on  the  repulsive  area  Arep  which  follows  the 
same  law  K,Niep  =  h'0E*Aiep  that  was  found  for  nonadhesive  contact  [4,  24]  with  a  similar 
value  of  k  ~  2.  Additionally,  the  attractive  load  is  given  by  the  average  force  y/Ac  times  the 
attractive  area,  where  Ac  is  the  range  at  which  the  potential  energy  drops  to  zero.  Hence 
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we  have 

JV  =  ^W«p  -  -t-Am  (32) 

K  A  c 

for  the  total  load  N. 

We  now  relate  the  attractive  area  Aatt  to  the  repulsive  area  Arep.  In  what  follows,  we 
assume  that  the  adhesive  interaction  only  weakly  perturbes  the  shape  of  the  gap  such  that 
we  can  use  the  3/2  power-law  separation  to  estimate  the  attractive  area.  Independent  of 
the  exact  functional  form  of  the  separation  of  the  contact  lips  g(R),  roughly  all  sites  within 
R  =  g_1(Ac)  around  the  repulsive  clusters  constitute  the  area  Aatt.  For  small  w,  large  Ac 
and  small  h'0  this  will  lead  to  mutually  overlaping  regions.  If  on  the  other  hand  w  3>  R  then 
the  attractive  region  will  be  small  compared  to  the  cluster  size  and  we  expect  Aatt  ~  PR- 
We  hence  estimate  the  attractive  area  to  be 


A 


att 


/  2AC 


A 


rep 


\h'0w 


(33) 


with  a  constantof  proportionality  f3  that  is  of  order  1.6.  Fig.  4  shows  the  value  of  f3  as 
obtained  from  the  numerical  calculations.  As  As  increases  the  value  of  f3  converges  towards 
roughly  unity.  For  small  As  we  expect  deviations  from  the  ideal  behavior  since  w  becomes 
of  order  of  the  A c/h'Q. 

The  value  of  (3  is  remarkably  indepedendent  of  the  interaction  strength  7.  Fig.  5a  shows 
the  dependence  of  this  factor  for  value  of  7  ranging  from  zero  to  0.1.  For  As  =  16  and 
Xs  =  32  this  value  is  remains  constant  at  roughly  1.6.  The  ratio  of  Aatt/Aatt,  however,  starts 
to  deviate  from  the  average  force  7/ A c  over  the  range  of  adhesive  energies  shown.  This  is 
most  likely  attributable  to  a  deviation  from  the  3/2  power  law  behavior  of  the  air  gap  due 
to  the  influence  of  the  adhesive  interaction. 

We  now  have  an  approximate  expression  for  the  total  load  load  N  as  a  function  of  the 
the  repulsive  (contact)  area  Arep.  Hence,  we  can  define  a  renomalized  k!  by 


I  =  — — —  =  -  — 4m  ( ?M s  ,34, 

k'  h'0E*Arep  k  h'0E* Ac  \h'0w ) 

that  depends  explicitly  on  the  adhesive  energy  7,  the  interaction  range  Ac  of  the  interatomic 
potential,  the  root  mean  square  slope  h'0  as  well  as  the  internal  length  scale  w.  Hence,  the 
influence  of  adhesion  is  controlled  by  the  short  wavelength  of  the  rough  surface  only. 

The  renormalized  k'  can  be  again  formally  split  into  a  repulsive  and  an  attractive  con¬ 
tribution.  Fig.  6  shows  the  dependence  of  these  individual  contributions  on  w.  Here,  the 
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FIG.  4:  (a)  Relationship  of  perimeter  P,  w  and  repulsive  area  vRep  (b)  Behavior  of  the  width  R  is 
the  attractive  area  as  a  function  of  w.  Here  R  is  the  average  R  over  the  whole  contact  geometry, 
(c)  Relationship  of  the  attractive  area  Hatt  to  perimeter  P  and  the  average  width  R.  (d)  Scaling 
of  Hatt/Hrep.  This  plot  is  the  product  of  (a),  (b)  and  (c).  In  all  plots,  symbols  denote  surfaces 
with  different  Hurst  exponents  H  and  root  nifj^n  square  slopes  h'0,  while  color  denotes  the  short 
wavelength  cutoff  As.  Symbol  size  is  a  measure  of  Al/As  where  only  data  with  Al/As  >  128  is 
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FIG.  5:  (a)  Dependence  of  the  ratio  of  attractive  to  repulsive  area  as  a  function  of  the  strength  of 
the  adhesive  interaction  7.  (b)  Dependence  of  the  ratio  of  attractive  load  to  attractive  area  on  7. 
For  larger  7,  this  ratio  starts  to  deviate  from  the  average  force  7/Ac. 


attractive  k  in  Fig.  6a  nicely  follows  the  relationship  Eq.  (34).  The  repulsive  k  (Fig.  6b)  is 
close  to  a  k  obtained  from  a  fully  nonadhesive  calculation  (Fig.  6c). 

Equation  (34)  now  allows  to  determine  the  onset  of  the  macroscopic  adhesive  behavior, 
i.e.  non-zero  pulloff  forces  for  the  two  surfaces.  Such  a  condition  is  obtained  when  N  <  0 
and  hence  1/k'  <  0.  From  this  condition  we  get 


w  <  2AC 


(«/?)* 

W>)f 


(35) 


as  the  upper  bound  on  the  admissible  values  of  w.  Since  w  ~  As  macroscopic  adhesive 
behavior  can  only  be  obtained  if  the  roughness  extends  down  the  length  scale  given  by 
Eq.  (35).  Equation  (35)  clarifies  the  result  of  Fuller  and  Tabor  [13]  that  roughness  can 
destroy  adhesive  interaction. 

From  this  scaling  relation  we  can  also  obtain  a  criterion  for  when  macroscopic  adhesive 
behavior  can  be  ruled  out.  The  value  w  =  1  is  a  natural  lower  bound  since  roughness  cannot 
extend  below  the  atomic  scale.  Then  we  find 


ao  (hp)  3  /  AA  3 
2§  k/3  \ao  J 


(36) 
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FIG.  6:  (a)  Attractive  k.  (b)  Repulsive  k.  (c)  Iv  obtained  from  a  nonadhesive  calculation  (i.e.  for 
7  =  0). 
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for  surfaces  to  become  adhesive.  All  dimensionless  quantities  in  Eq.  (36)  are  on  the  order  of 
unity,  giving  roughly  5/a0  >  (h'0)  I . 


[1]  H.  Hertz,  J.  reine  angew.  Math.  92,  156  (1881). 

[2]  K.  L.  Johnson,  Contact  Mechanics  (Cambridge  University  Press,  1985). 

[3]  B.  Luan  and  M.  O.  Robbins,  Nature  435,  929  (2005). 

[4]  S.  Akarapu,  T.  Sharp,  and  M.  O.  Robbins,  Phys.  Rev.  Lett.  106,  204301  (2011). 

[5]  L.  M.  Falicov  and  F.  Yndurain,  J.  Phys.  C:  Solid  State  Phys.  8,  147  (1975). 

[6]  V.  R.  Velasco  and  F.  Yndurain,  Surf.  Sci.  85,  107  (1979). 

[7]  Y.  Saito,  J.  Phys.  Soc.  Jpn.  73,  1816  (2004). 

[8]  C.  Campana  and  M.  H.  Miiser,  Phys.  Rev.  B  74,  075420  (2006). 

[9]  C.  Creton  and  S.  N.  Gorb,  MRS  Bulletin  32,  466  (2007). 

[10]  K.  Autumn,  M.  Sitti,  Y.  A.  Liang,  A.  M.  Peattie,  W.  R.  Hansen,  S.  Sponberg,  T.  W.  Kenny, 
R.  Fearing,  J.  N.  Israelachvili,  and  R.  J.  Full,  Proc.  Natl.  Acad.  Sci.  USA  99,  12252  (2002). 

[11]  K.  Autumn,  Y.  A.  Liang,  S.  T.  Hsieh,  W.  Zesch,  W.  P.  Chan,  T.  W.  Kenny,  R.  Fearing,  and 
R.  J.  Full,  Nature  405,  681  (2000). 

[12]  P.  F.  A.  Maderson,  Nature  203,  780  (1964). 

[13]  K.  N.  G.  Fuller  and  D.  Tabor,  Proc.  R.  Soc.  Lond.  A  345,  327  (1975). 

[14]  M.  Lane,  Annual  Review  of  Materials  Research  33,  29  (2003). 

[15]  J.  N.  Israelachvili,  Intermolecular  and  Surface  Forces  (Academic  Press,  London,  1991). 

[16]  G.  Huber,  H.  Mantz,  R.  Spolenak,  K.  Mecke,  K.  Jacobs,  S.  N.  Gorb,  and  E.  Arzt,  Proc.  Natl. 
Acad.  Sci.  USA  102,  16293  (2005). 

[17]  A.  R.  C.  Baljon  and  M.  O.  Robbins,  Science  271,  482  (1996). 

[18]  C.  Gay  and  L.  Leibler,  Phys.  Rev.  Lett.  82,  936  (1999). 

[19]  D.  S.  Dugdale,  J.  Mech.  Phys.  Solids  8,  100  (1960). 

[20]  D.  Maugis,  J.  Colloid  Interf.  Sci.  150,  243  (1992). 

[21]  B.  Gotsmann  and  M.  A.  Lantz,  Phys.  Rev.  Lett.  101,  125501  (2008). 

[22]  B.  Luan  and  M.  O.  Robbins,  Phys.  Rev.  E  74,  26111  (2006). 

[23]  B.  B.  Mandelbrot,  The  Fractal  Geometry  of  Nature  (W.  H.  Freeman,  1982). 

[24]  S.  Hyun,  L.  Pei,  J.  Molinari,  and  M.  O.  Robbins,  Phys.  Rev.  E  70,  26117  (2004). 


18 


